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ABSTRACT 


The evolution of the spin rate of comet 9P/Tempel 1 through two perihelion 
passages (CYs 2000 and 2005) is determined from 1922 Earth-based 
observations taken over a period of 13y as part of a World-Wide observing 
campaign and 2888 observations taken over a period of 50d from the Deep 
Impact spacecraft. We determine the following sidereal spin rates (periods): 
209.023± 0.025 7day (41 .335 ± 0.005 h) prior to the 2000 perihelion passage, 
210.448±0.016 7day (41.055 ± 0.003 h) for the interval between the 2000 and 
2005 perihelion passages, 21 1 ,856± 0.030 7day (40.783± 0.006 h) from Deep 
Impact photometry just prior to the 2005 perihelion passage, and 21 1 ,625± 0.012 
7day (40.827± 0.002 h) in the interval 2006-2010 following the 2005 perihelion 
passage. The period decreased by 16.8 ± 0.3 min during the 2000 passage and 
by 13.7 ± 0.2 min during the 2005 passage suggesting a secular decrease in the 
net torque. The change in spin rate is asymmetric with respect to perihelion with 
the maximum net torque being applied on approach to perihelion. The Deep 
Impact data alone show that the spin rate was increasing at a rate of 
0.024±0.003 7d/d at JD2453530. 60510 ( i.e 25.134 d before impact) and 
provides independent confirmation of the change seen in the Earth-based 
observations. 

The rotational phase of the nucleus at times before and after each 
perihelion and at the Deep Impact encounter is estimated based on the Thomas 
et al. pole and longitude system (2007, Icarus 187, 4-15). The possibility of a 
180° error in the rotational phase is assessed and found to be significant. 
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Analytical and physical modeling of the behavior of the spin rate through of each 
perihelion is presented and used as a basis to predict the rotational state of the 
nucleus at the time of the nominal (i.e., prior to February 2010) Stardust-NExT 
encounter on 201 1 February 14 20:42. 

We find that a net torque in the range of 0.3 -2.5x1 0 7 kg.m 2 .s' 2 acts on 
the nucleus during perihelion passage. The spin rate initially slows down on 
approach to perihelion and then passes through a minimum. It then accelerates 
rapidly as it passes through perihelion eventually reaching a maximum post- 
perihelion. It then decreases to a stable value as the nucleus moves away from 
the sun. We find that the pole direction is unlikely to precess by more than 
~1°/perihelion passage. The trend of the period with time and the fact that the 
modeled peak torque that occurs before perihelion is in agreement with 
published accounts of trends in water production rate and suggests that 
widespread H 2 O out-gassing from the surface is largely responsible for the 
observed spin-up. 
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1. Introduction. 


In their assessment of the spin state of 9P/Tempel 1 following the Deep Impact 
encounter A’Hearn et at. (2005) found a sidereal period of 1 .701 ± 0.014 d (21 1 .6 
±1 .7 7d). Thomas et al., (2007) determined the sense of spin to be direct and the 
preliminary pole position was RA = 5°, Dec = +78° (± 10° on the sky) roughly 1 1° 
from one of the two, photometrically degenerate, directions (46°, +73°) found 
earlier by Belton et al. (2005) using ground-based photometry. A’Hearn et al. also 
noted that the pre-impact rotation period of 1 .744 ± 0.006 days determined by 
Belton et al. (2005) differed significantly from the spacecraft value and suggested 
that the difference might be explained by an inadvertent shift in the analysis of 
earth-based data by a half or whole cycle between observing runs. Later, Belton 
et al. (2006) improved the Deep Impact sidereal period estimate to 1 .6976 ± 
0.0096 d (21 2.06 ± 1 .2 °/d) and Thomas et al (2007) revised the pole position 
(J2000) to RA = 294°, Dec = 73° (± 5° on the sky). Thomas et al. also set the 
prime meridian as W(t) = 252.63° + 212.064° t, where t is the number of days 
since the standard epoch (JD 2451545.0). W(t) is the angle between the chosen 
prime meridian and the intersection of the body equator and the standard Earth 
equator and defines the rotational phase of the nucleus at time t. This latter 
formula assumes a constant rotation period of 1 .6976 d between the time of 
impact (JD 2453555.73928) and the standard epoch. 

The small difference between the pre-impact rotation rate and the 
spacecraft value would have been of little concern had it not been for the 
recognition by J. Veverka and his colleagues that the Stardust spacecraft, which 
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had recently encountered comet 81 P/Wild 2 (Brownlee et al., 2004), was 
hibernating in deep space, could be revived, and had enough propulsion 
capability to reach 9P/Tempel 1 for an encounter on 201 1 February 14. This 
mission, now called Stardust-NExT, was selected as a Discovery mission of 
opportunity by NASA ( www.astro.cornell.edu/next/Science.htm ). A Level 1 
science requirement of this mission is to “Image 25% of the surface previously 
observed in the Deep Impact mission at better than 80 m/pixel” in order to look 
for changes in the condition of the surface that might have occurred during the 
previous perihelion passage (201 1 January 12.2). A secondary science goal is to 
image the, as yet unseen, artificial crater formed during the Deep Impact mission 
that is located at 350. 4W, -29.1 (Thomas, 2010, private communication). To 
ensure that these objectives can be met, a high-precision rotational ephemeris 
and an assessment of its stability is required and it is for this reason that the 
present study was initiated. 

In subsequent preparations for the Stardust-NExT mission it was noted 
that the Deep Impact rotation rate calculated by Belton et al. (2006) did not 
correctly phase the light curves obtained some 14 months earlier from the 
Hubble and Spitzer Space telescopes (Lamy et al., 2007; Lisse et al., 2005). 

This was the first quantitative indication that comet 9P’s rotation might be 
changing as it approached perihelion. 

The theoretical basis for short timescale changes in cometary spin has 
been emphasized by Jewitt (1997; also summarized in Jewitt, 2004; see also 
Samarasinha et al. 2004) and exploratory calculations of excitation timescales 
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have been carried out by Gutierrez et al. (2002), Jorda and Gutierrez (2002), and 
Gutierrez and Davidsson (2007). For a small (effective radius = 3.0 ± 0.1 km; 
Thomas et al., 2007), underdense (bulk density ~ 400 kg.m' 3 ; Richardson et al. 
2007) nucleus with a water production rate of 6 x 10 27 molecules/s (Schleicher et 
al., 2006) the timescale for substantial changes in the spin state is ~ 90 y based 
on Jewitt’s formulation of spin-up time and his conjecture that the typical 
dimensionless moment arm for torques is -0.05. Thus, from a theoretical point of 
view it should not be surprising if 9P/Tempel 1 was changing its current period by 
- 1% (0.4 h) in a single perihelion pass or if the direction of the rotation pole 
drifted by a degree or two. 

There is also a growing observational base to support the measurable 
presence of this effect in comets. Drahus and Waniak (2006) have shown 
through the introduction of a novel photometric time-series analysis technique 
that the rotation rate of the distant comet C/2001 K5 (LINEAR) was perceptibly 
spinning-down as it receded from perihelion passage. In addition, earlier studies 
have found evidence of possible changes in spin rate in comets lOP/Tempel 2 
and 6P/d’Arrest (Mueller and Ferrin, 1996; Gutierrez et al., 2003). Other evidence 
of the action of rotational torques includes the cases of comets 1 P/Halley, 
2P/Encke and 29P/Schwassmann-Wachmann 1 each of which have been found 
to be in rotationally excited states (Belton et al., 1991; Samarasinha and 
A’Hearn, 1991; Meech et al., 1993; Belton et al., 2005). 

Deep Impact photometry and imaging data from the ongoing Worldwide 
Earth-based campaign on 9P/Tempel 1 (Meech et al, 2005, 2011) plus an early 
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Hubble Space Telescope study by Lamy et al. (2001) provide an unprecedented 
set of data with which to investigate the stability of the spin state of 9P/Tempel 1 . 
The data that we use from ground-based and HST sources are described in 
Section 2 where we separate them into three groups: Region A (1997-1999), 
Region B (2001-2004) and Region C (2006 - 2010). This allows us to document 
the changes that occurred during the 2000 and 2005 perihelion passages. In 
Section 3 we present the Deep Impact approach photometry that we use to 
obtain direct evidence for an acceleration of the spin rate. In Section 4 we 
provide the theoretical basis and assumptions used in the analysis of the data. In 
section 5 we outline the rotational analysis and present the basic results on the 
spin rate of the nucleus and its rotational phase. In section 6 we discuss the 
dynamical evolution of the comet’s spin state and construct analytical models for 
its changes through perihelion passage. In Section 7, we provide a general 
discussion of the relationship of our results with previously published studies of 
the comet’s H 2 0 production rate. We also use our results to predict the rotation 
state we expect will be experienced by the Stardust-NExT mission at its 
encounter with 9P/Tempel 1 on14 Feb, 201 1 . Section 8 contains a summary of 
our primary conclusions. 


10 



2. Earth-based observations and two independent methods of analysis. 

In an accompanying paper, Meech et al. (2011) provide a detailed description of 
the Deep Impact World-Wide campaign, its goals, participant contributions, 
observations and results. In Fig. 1 we show R(1 ,1 ,a) magnitudes for the entire 
data set after reduction to unit heliocentric and geocentric distance and where a 
is the solar phase angle. The rise and fall of the coma brightness around 
perihelion dominates the figure and the substantial effect of diminishing solar 
phase angle at the oppositions of 2000, 2002, 2003, 2004, 2006, 2008 and 2010 
is clear. A small sub-set of this data was obtained as V magnitudes and these 
have been converted to R magnitudes using (V-R) = 0.50 mag (Li et al., 2007). 
Embedded in this data set are the results of three studies done with the Hubble 
Space Telescope ( HST ). These were obtained in 1997 (Lamy et al., 2001), 2004 
(Lamy et al., 2007) and 2009 (Meech et al., 2011). The latter two were done at 
the request of the Deep Impact team in order to obtain data of sufficient quality to 
distinguish between alias periodicities and to measure, with the highest possible 
accuracy, the rotational phase of the nucleus at each epoch. The V(1,1,a) 
magnitudes in Lamy et al. (2007) were converted to R(1 ,1 ,a) as noted above. 

All of the above photometric reductions depend on orbital information 
obtained from the JPL Horizons site ( http://ssd.ipl.nasa.gov/7horizons ) using the 
default orbit solution (K054/15) for 9P/Tempel 1 . The complete set of reduced 
ground-based data, relevant geometry, and the timing used in this paper is listed 
in a supplementary electronic data file appended to Meech et al. (201 1 ). In our 
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analysis we divided the data into three intervals each of which was expected to 
yield significantly different values for the spin rate: Region A, before the year 
2000 perihelion passage; Region B, between the 2000 and 2005 perihelion 
passages; and Region C, post the 2005 perihelion passage. A fourth region D 
was reserved to cover the time post the 201 1 perihelion passage and through the 
Stardust-NExT encounter on 14 Feb, 2011. 

Correction for solar phase angle brightness effects. Li et al. (2007) have 
determined the disk integrated phase function of 9P to be p = 0.046 ± 0.007 
mag./deg for 4° < a < 1 1 7°. Earlier Belton et al. (2005), using data obtained 
between 1997 and 2002, found evidence for an increase in (3 inside of 4°. They 
represented the phase law as a polynomial (AR(1,1,a) = -0.0180955 - 0.250260a 
+ 0.0306201a 2 - 0021 805a 3 + 0.0000798a 4 - 0.0000015a 5 mag.) good for a < 15°. 
In Fig. 2 the data and these two phase laws are compared, but, as can be seen, 
a simple linear regression, R(1,1,a) = 14.905 + 0.0449a mag, gives an excellent 
account of the current data for 1° < a < 14°. It is this latter relationship that we 
have used to correct the magnitude data to zero solar phase angle (Fig. 3). 

In Figs. 3, 4 and 5 we show R(1 ,1,0) as a function of calendar year, time 
from perihelion passage, and heliocentric distance. In Fig. 3 we see how well the 
solar phase angle brightness effect has been removed. The scatter in 
magnitudes while the comet is near aphelion now primarily reflects the variation 
in brightness caused by the spin of 9P’s irregularly shaped nucleus. In Fig. 4 we 
see how consistent the data are from one perihelion passage to the next. In Fig. 
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5 we see how the brightness of the coma is maintained at higher levels post- 
perihelion relative to its behavior pre-perihelion, an effect seen in other comet 
light curves and that prevented our use of early post-perihelion observations in 
the determination of the rotational state. 

Methods and preparation of the data for rotational analysis. Because of our 
intention to use the results of this study to adjust the arrival time of the Stardust- 
NExT spacecraft at the comet (and so meet Level 1 requirements specified by 
NASA) it was decided to perform two independent analyses for the predicted 
rotational state of the comet. While the details of each of these are reserved to a 
later section, we now give an overview of these independent techniques as an 
introduction. 

The first was done at the Jet Propulsion Laboratory in Pasadena and was 
based on a least-squares fit of model light curves to each region of R(1 ,1 ,0). The 
model light curves were generated using a combination of the Thomas et al. 
(2007) shape model, a Hapke photometric function, SPICE orbital data 
( http://naif.ipl.nasa.gov/naif/data comet.html ), and an unexcited rotation model, 
i.e., the spin axis was assumed fixed in space and coincident with the principal 
axis of maximum moment of inertia. Hapke parameters for 9P have been 
determined by Li et al. (2007), however, in certain applications that utilized 
intensive computations a Lommel-Seeliger (LS) function was substituted. The LS 
photometric function represents single (isotropic) scattering from a semi-infinite 
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medium and comparisons with Hapke calculations demonstrated that it was an 
excellent proxy for 9P which has a low surface albedo. 

The second investigation, which was done in parallel with that at JPL, was 
done at Belton Space Exploration Initiatives, LLC, in Tucson. This study applied 
standard astronomical period finding techniques to R(1 , 1 ,0) after removing orbital 
synodic and solar phase angle timing (Harris et al., 1984) effects. The data set 
was transformed to one that would have been acquired by an inertial observer 
fixed relative to the comet. This analysis avoided the use of the shape model and 
the choice of a surface photometric function. Under normal circumstances 
cometary light curve data rarely extend over an interval longer than two or three 
months near a single opposition and such observations can be phased to 
determine an adequate approximation to the rotational period without first 
accounting for orbital synodic effects or worrying about a change in the timing of 
light curve maxima due to changing illumination geometry (Fig. 6). In fact, these 
corrections would be hard to make without prior knowledge of some of the 
physical properties of the nucleus. In the present case, such prior knowledge is 
available, i.e., the rotational period, the polar axis, the shape, and photometric 
behavior of the nucleus are all approximately known (Belton et al., 2005; Thomas 
et al., 2007; Li et al., 2007) and first order corrections can be made. Since the 
observations are, as shown in Fig. 7, spread over several oppositions that are 
widely separated around the orbit and involve a wide range of solar phase 
angles, accounting for these effects is essential to obtaining well-defined light 
curves and in achieving the highest accuracy in the sidereal spin rate. 
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Both of these investigations were subjected to periodic and independent 
peer reviews during the analysis. The internal peer review team included A. 
Harris, T. Duxbury and D. Scheeres. 

3. Deep Impact approach photometry. 

As noted in the introduction a number of estimates have been made of 
the rotational period of the nucleus based on the Deep Impact approach 
photometry. However, all of these where based on an early form of the 
photometry that was subsequently found to have short-comings that may have 
affected the results in the earliest parts of the approach sequence. These 
problems, which could possibly affect the accuracy of the estimated period, 
include estimation of the bias correction to the nearest DN (data number, a linear 
measure of the brightness), ignoring faint horizontal striping in the images, and 
ignoring small corrections to the photometry that are required when the comet 
happened to be placed on the two rows in the image surrounding the horizontal 
boundary between the upper and lower halves of the detector. These corrections 
have now been made (see appendix A for a more complete discussion) and a 
new, improved, version of the approach photometry made available for analysis. 

It consists of 595 points of “science” data covering 63. Id on approach and 2419 
points of “navigation” data covering an interval of 49.8 d. The last datum in these 
sets was taken at 0.25 d before impact. All data considered here (2888 Science 
and Navigation observations taken from 50 to 0.25 d before impact) were taken 
with the Medium Resolution Instrument (MRI) through one or the other of two 
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CLEAR, effectively identical, filters with a central wavelength at 650 nm 
(Hampton et al., 2005). The photometry, which refers to circular apertures 5, 7, 

9, 15, 20, 25, and 30 pixels in diameter centered on the comet, can be found on 
the PDS Small bodies Node http://pdssbn.astro.umd.edu in the Deep Impact 
archives. 

We assume that the light of nucleus, which was centered in the circular 
aperture, was completely contained within it. At 0.25 d before impact (the last 
datum) the mean diameter of the nucleus subtends 1 .8 pixels; i.e., the nucleus is 
always within our smallest aperture of 5 pixels. The DN at each time, t, consists 
of two parts: that contributed by the nucleus and that by the inner coma. We 
express the signal as: 

DN (t) = p.Cn. fn(a). F n (t). r - 2 .d ' 2 + F c (t, r, d, a) 

Where (3 is a constant calibration factor, C n is the mean brightness of the 
nucleus, f n is the solar phase function (as determined in Section 2). F n (t) is the 
variation of the brightness of the nucleus as it rotates and r and d are the 
distances from the sun and the spacecraft in AU. F c is the coma contribution to 
the signal and has a complex dependence on the parameters shown. With the 
exception of p, F c and F n all of the other quantities are specified. Our objective is 
to determine the product p.F n and analyze its time dependency for the rotational 
period. To do this we must first separate p.F n from F c , which we assume to make 
a negligible rotational contribution to the signal. Unfortunately, this latter 
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assumption cannot be precisely true since the origin of coma material is tied to 
the surface of the nucleus and will initially share its motion. As the coma material 
flows out from the nucleus conservation of angular momentum will reduce its 
angular motion and periodicities may be introduced that are systematically 
different to that of the nucleus. 

We have investigated the possible effect of this phenomenon by searching 
for coma variability in the difference signal between the 7 and 5 pixel apertures. 
This signal should arise entirely from the coma. For times close to impact, when 
the very inner coma is being sampled, we have detected variability at the nucleus 
period but with an amplitude that is -0.1 that of the variability of the nucleus itself. 
At 10 days before impact, when a more distant and broader region of the coma is 
being sampled, coma variability is undetectable presumably being overwhelmed 
by noise. At 0.25 d before impact the coma signal in the 5 pixel aperture is only ~ 
0.04 that of the nucleus (Fig. 8) and so the amplitude of the coma’s observed 
variability contributes a negligible ~ 0.004 fraction of the total observed 
amplitude. At earlier times, well before impact and when the coma is the 
dominant component of the signal, the 5 pixel aperture integrates over a wide 
region of the coma presumably washing out any rotational variability from that 
source. In Fig. 8 we show the approximate contribution of the nucleus to the total 
signal in the 5 pixel diameter aperture. This first order separation of the nucleus 
from the coma was done by assuming that there is a linear relationship between 
the coma signal and aperture diameter in apertures 5 through 30. The signal 
extrapolated to D = 0 is a first order estimate of the contribution of the nucleus. 


17 



This simple coma model is based on a symmetric constant velocity outflow in the 
coma and a point nucleus and predicts the roughly linear dependence of the 
signal with aperture size that is seen in the data. The coma is seen to dominate 
until about 4 days before encounter and considerable care was taken to remove 
it before the rotational analysis. 

The data consists of two groups: “Science” data and “Navigation” data 
whose time coverage, sampling, and mode of acquisition were quite distinct. 
While they are taken with the same instrument and filters they were taken in 
different camera modes as described in appendix A. Nevertheless, the two types 
of data were found to be photometrically consistent and we use them together as 
a “joint” data set. To illustrate the separation of coma from (3.F n and the 
“cleaned-up” data we plot, in Fig. 9, the normalized brightness, r 2 .d 2 .DN/C n .f(a), 
versus time from impact. C n is taken as 100 units and r, d and a were obtained 
from the JPL Horizons system where the Deep Impact spacecraft is identified as 
“@ - 140”. The coma and mean nucleus is represented as a 4 th -order polynomial 
in each of the two cleaned-up data sets and then subtracted leaving the time 
variable component of the nucleus, (3.F n . In principle the mean level of this 
component should be constant and zero and this can be seen to be 
approximately the case in Fig. 9. The two data sets are then combined to form 
the joint version of p.F n shown in the lower panel of the figure. 


4. Rotational equations of motion and assumptions 
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The action of forces generated by the momentum of gas and dust leaving the 
surface of an active comet nucleus can at any instant of time be decomposed 
into two parts: those which act at the center of mass of the body, F, and those 
which apply torques, T. Both forces are functions of time, t, and vary rapidly on 
a rotational time scale but are thought to vary relatively smoothly on an orbital 
time scale. 

Occasionally major cometary events such as a splitting, or a major 
outburst, or the appearance or disappearance of a major active region, or a close 
encounter with a major solar system body may occur that lead to unpredictable 
changes in these forces. Such effects are not considered here and we assume 
that the rate of change of the average of T(t) over a rotational cycle, d<T(t)>/dt, 
changes smoothly over orbital timescales. It is not necessary that such orbital 
changes are the same from one orbit to another and we shall, in fact, find that 
they are not. 

We assume that the nucleus rotates as a rigid body and the vector 
equation for the angular motion of such a body in a fixed frame, e.g., Rutherford, 
(1951), is: 

dhldt = T tidal (t) + T(t) (1) 

Where h is the angular momentum of the body and T ti dai are torques that come 
into play during near encounters with planets or the sun. In the present 
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application we shall assume that T t idai = 0. The net torque due to mass loss is 
therefore defined as: 

T = J(p x q.V).dS (2) 

where p is the position vector of an elemental area dS with respect to the center 
of mass of the nucleus, V(t) is the velocity of the outflow and q(t) is the net rate 
of mass loss from dS. The integration is over the entire surface and, although not 
specifically indicated in the notation, averaged over a rotational cycle. With co as 
the angular velocity of the nucleus relative to fixed axes that momentarily 
coincide with the principle axes of inertia in the nucleus, Eq.1 yields Euler’s 
equations of angular motion: 

A*dco x /df + (C*-B*) . COyCOz = [f(p x q.V).dS] x 

B'dayldt + (A*-C*).ro z (Ox = [J(p x q.VJ.dS], (3) 

C*dco z /df + (B*-A*) . COxCOy = [j(pxq.V).dS] z 

Where A*, B*, and C* are the principal moments of inertia. 

The assumption that a cometary nucleus rotates as a rigid body deserved 
some comment in light of the low cohesive strength found for cometary material 
in D/Shoemaker-Levy 9 (Asphaug and Benz, 1996) and 9P (Richardson eta!., 
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2007) and the proposition of Belton and Melosh (2009) that fluidized transport of 
substantial amounts of material may take place episodically in the interior 
ultimately leading to outbursts and surface flows. To justify our assumption it is 
sufficient to show that rotation is a minor source of stress and that the loss of 
material in a repetitive outburst has a minor effect on the moment of inertia. Near 
the ends of the long axis of 9P the gravity is -0.027 cm/s 2 and the centripetal 
force is - -0.001 cm/s 2 (Thomas, 2007; private communication) thus the rotational 
contribution is minor. We estimate the moment of inertia of the nucleus to lie 
between 0.8 - 4.6 x 10 19 kg.m 2 assuming a homogeneous mass distribution in 
the interior, a spherical approximation to the shape, and using the mass range 
determined by Richardson et al. (2007). A typical repetitive outburst releases 
about 10 6 kg of material at the surface whose maximum contribution to the 
moment of inertia for a 3 km radius body is - 10 13 kg.m 2 . Since this is much less 
than the moment of inertia, the loss of material in a repetitive outburst, or similar 
events, will not have a noticeable effect on the rotational dynamics of the 
nucleus. 

These equations provide a firm basis for detailed ab initio simulations of 
the evolution of the spin of the nucleus, if the inertias and mass loss are fully 
understood e.g., Samarasinha and Belton (1995) or Gutierrez and Davidsson 
(2007). The study of non-gravitational effects on 9P’s orbit suggest that the 
direction of the spin axis is stable or, at most, slowly changing (Yeomans et al., 
2004b). In addition, the close similarity of the pole direction (RA, Dec = 294°, 
+73°; J2000) found by Thomas et al. (2007) from Deep Impact images of the 
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resolved nucleus with the pole direction (RA, Dec = 293°, =+73°) derived from 
the time dependence of coma features by Vincent et al. (2010) suggest that the 
direction of the polar axis is well determined. In addition the pole direction (RA, 
Dec = 46°, +73 revised to 317°, +81° in 2006) derived from light curves obtained 
over several years prior to the Deep Impact encounter (Belton et al., 2005) 
support the idea that any precessional drift in the pole direction is slow, i.e., not 
exceeding a few degrees on the sky per perihelion passage. Moreover, with the 
assumption of a homogeneous mass distribution in the nucleus, observations of 
the shape of the nucleus allow us to estimate the direction of the principal axis of 
maximum moment of inertia, which, within the errors of estimation, is found to be 
parallel and coincident with the estimated spin axis (Thomas et al. 2007). All of 
these observations support the idea that the spin of the comet is close to its fully 
relaxed state. Additional support comes from an analysis of periodicities in the 
Deep Impact light curve, which is found to yield only harmonics of a single 
period. While this does not ensure that the nucleus is in a fully relaxed state it is, 
nevertheless, a necessary condition. We conclude that the nucleus is apparently 
close to the state of simple rotation around its principal axis of maximum moment 
of inertia, which we take as C*. In this case the z-axis is coincident with the spin 
axis and we assume that co x = co y = dco x /dt = dco/dt = 0. This assumption implies 
that, for over the period of the observations being considered, the direction of the 
rotation axis is considered fixed. We will revisit and provide a check on this 
assumption in Section 6. 

After dropping the z subscript, Eqs. 3 reduce to: 
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cfa>(t)/df= T(t) 


( 4 ) 


where 


T(t) = J(p x q.V).dS /C* 
averaged over a rotational cycle. 

We define the rotational phase, W(t), of the nucleus at time, t, as the 
angular distance of the prime meridian, as defined by Thomas et al. (2007), to 
the meridian that contains the direction of the ascending node (RA, Dec = 23.80°, 
0.0°; J2000) as seen from the nucleus, i.e., 

dW(t)/dt = co(t) (5) 

Integrating from t 0 to t we get 


W(t) = W 0 + lto,tco(t).dt 


( 6 ) 


Since the sense of spin of 9P is observed to be direct (Thomas et al., 2007), W(t) 
is the West longitude of the reference direction at time t. At the time of impact in 
the Deep Impact Mission (JD2453555. 73928) Thomas et al. found that 

W(t) = 252.63° + 21 2.064*(t- 1 0 ) (7) 
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where the standard epoch is t 0 = JD2451 545.0 and the assumed spin rate 
(212.064°/d) was constant between the standard epoch and the time of impact. 

In the analysis that follows we make two further assumptions. First, we 
assume that for an extended period around aphelion (Regions A, B, C, D) doVdt 
=0, i.e., torques during that period are negligible. Secondly, we assume that the 
light curve of the nucleus, in the same extended period around perihelion, can be 
predicted from Thomas et al’s. (2007) shape model. Referring to Eq. 6 and Fig. 

3, the first of these assumptions implies that the following relationships exist: 

Wj(t) = W 0j + Sj*(t - 1 0 ) j = Regions A, B, C, D (8) 

and the rotational analysis in the next section is designed to discover W 0 j and Sj 
for each of the regions. 

The second of the above assumptions is more problematical because the 
3-dimensional shape of 9P’s nucleus is poorly defined over large areas of the 
nucleus. This is a direct result of the linear Deep Impact flyby geometry and the 
slow rotation of 9P. Tests of model predictions with Deep Impact approach 
photometry show that with either a Lommel-Seeliger or Hapke scattering function 
the Thomas shape model can give a good match to the observations. However, 
while the shape model gives good results for the light curve at most rotational 
phases, we find that there is a restricted range of rotational phase (see section 5 
below) where the data and models are discordant. 


24 



The validity of both of the above assumptions depends on the premise 
that there is an extended interval near aphelion when the effect of coma on the 
light curve and the dynamics of the nucleus is negligible. To assess the influence 
of the coma we first remove the mean brightness of the nucleus from the data 
shown in Fig. 3. The mean absolute R magnitude of the nucleus at zero solar 
phase angle is taken as 14.905 (see Section 2 for the origin of this value) and 
magnitudes are converted to relative brightness units with the mean brightness of 
the nucleus set at 100 units. In Fig. 10 the pre- and post-2005 perihelion coma 
brightness is compared to the mean brightness of the nucleus (shown as a 
dashed horizontal line). Pre-perihelion the signal is essentially coma free at 
distances beyond 3.2 AU where the RMS variability of ±15 units is mainly due to 
rotation. The error of an observation is typically ± 3.5 units. In the post-perihelion 
period only the data beyond 4.1 AU can be considered free of coma. The RMS 
spread at these distances is again ±15 units. The post-perihelion data between 
3.6 and 3.9 AU is from the 2006 opposition and is ~ 42% due to coma. However, 
the RMS variation of ±17 units is roughly similar to the mean level found beyond 
4.1 AU suggesting that the variability in 2006 remains dominated by rotation. A 
substantial observational effort was made during the September 2000 opposition 
to obtain rotational information when the comet was near 3 au post-perihelion. 
However, at that time the coma level was 2.7 times that of the nucleus and the 
RMS spread had increased to ± 55 units. Since the signature of nucleus rotation 
is expected to be at the ± 15 unit level, it is clearly overwhelmed by variability in 
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the coma. We have therefore omitted this part of the data from the rotational 
analysis. 

5. Rotational analysis 

We determine the values of Woj and Sj (Eq. 8) for each region j using the two 
methods described earlier. They share identical data sets and the combined 
ground-based and HST data for each region are shown in Fig. 1 1 . The Deep 
Impact data are already displayed in the bottom panel of Fig. 9. Descriptive 
information on the data is collected in Tablel . 

5. 1 The JPL Method. This is simply a classical least squares fit of a rotation 
model to the measured photometry. For the more typical case of unaccelerated 
rotation, which is presumed between active periods, the model parameters are 
the rotation phase Woj=Wj(to) and rate Sj=Sj(to) at some epoch to. In this way the 
spin phase and rate were obtained from a simultaneous fit to the photometry, in 
contrast to the power spectrum analysis used in the Tucson method described 
below. 

The least squares approach attempts to minimize the sum of squares of 
the photometric residuals, which are the differences between the observed and 
computed (O-C) magnitudes of the cometary nucleus at each measurement time 
The observed photometry is reduced to R(1 , 1 ,0) values, while the computed 
magnitudes are obtained from a special-purpose synthetic light curve generation 
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tools (“runlcrv.cgi” and “runlcrv_hapke.cgi” provided by B. Carcich) that can be 
queried online. 

The synthetic light curve tool takes as input the rotation history of the 
comet based on the assumed values of Woj and Sj and returns the received flux 
at a list of requested observation times for a specified observing location. The 
position of the observer can be selected as the geocenter, Deep Impact 
spacecraft, Stardust-NExT spacecraft, Spitzer Space Telescope or Hubble 
Space Telescope. These tools use Lommel -Seeliger and Hapke scattering laws 
and the Thomas et al. (2007) shape model. They use a plate model for the comet 
shape to capture the illumination effects for the particular positions of the 
observer and the sun at the emit time, and the light time delay is fully 
incorporated in order to accurately represent the measured flux at the 
observation time, i.e., receive time. We used a plate model having vertices on 6° 
centers, which we found to be a suitable compromise between performance and 
fidelity. Fluxes were converted to generic magnitudes and the optimal offset to R- 
band was also estimated. Some individual batches of photometry also required 
the estimation of an independent, ad hoc magnitude offset, typically due to coma 
contamination or photometric calibration issues. 

The least squares fitting process converged well when near a local 
minimum, but nonlinearities often led to large corrections when the initial guess 
was far from the minimum. Thus special care was needed in slowly building up to 
a global fit to larger data sets, e.g., an entire quiescent period. A typical fitting 
approach was to estimate the spin state with a short (few to several days) but 
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relatively dense data set, where the light curve variation was clearly visible. The 
epoch to would be situated within this data set (Fig. 12). Adding more and more 
data to slowly extend the fit span occasionally yielded good results, but in many 
cases the measurement gaps were too large to prevent ambiguities in the 
number of intervening rotations. An alternate method was to fix the rotation 
phase obtained from the short, dense data set and then scan a wide range of 
rotation rates to see which ones fit the larger data set the best. This approach 
yielded something akin to the periodograms used in the Tucson approach, but 
which had a somewhat different origin and interpretation (Fig. 13). With this 
approach we could generally identify a single or perhaps a few candidate 
frequencies that could be considered more carefully. 

Overall, fits were generally acceptable, but as can be seen from Fig. 17, 
the model light curve did have a significant departure from the measured light 
curve. This was associated with illumination of regions where the shape model is 
not well constrained by the Deep Impact approach photometry. However, the 
fitting process need only fit the gross features of the light curve and so 
mismodeling of the light curve was not a serious obstacle to estimating the 
rotation state, aside from the 180° phase ambiguity discussed in Sec. 5.2 below. 

Ground-based Photometry. This large data set, described more 
exhaustively by Meech et al. (201 1 ), formed the foundation for all of the 
quiescent fits in Regions A, B, and C. Most of the photometry obtained while the 
comet was highly active was not usable for rotation estimation, but many good 
batches even with substantial coma contamination were found to be helpful. For 
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the JPL method, isolated individual observations contributed little to rotation 
knowledge, but the numerous batches ( e.g Fig. 12) showing clear light curve 
variation and slope from night tonight proved to be vital. 

Deep Impact Photometry. The densest set of photometry available for this 
study was that derived from the approach photometry from the Deep Impact 
mission. This data set followed the model light curve with excellent coverage 
over seven weeks (Fig. 14) and, moreover, it was obtained while the comet was 
active and the spin state was presumably accelerating. With this in mind, we 
extended our rotation model to estimate the angular acceleration of the comet 
during the Deep Impact approach and found that the comet was indeed spinning 
up during this time. The acceleration model applied a torque proportional to the 
sublimation rate of water, as given by the g(r) function commonly used to model 
nongravitational accelerations on comets (Marsden et al. 1973). This approach 
allows the acceleration to build slowly, reaching a peak at perihelion and then 
fading back to effectively zero around 2.5 AU post-perihelion. As we explain 
beow, the comet acceleration profile was rather more complex, but we were still 
able to use this simpler model to estimate the acceleration on the relatively short 
interval of the Deep Impact approach. Specifically, we find that the best-fitting 
acceleration according to this model is dS/dt = A * g(r), where A = 0.0789±0.0031 
deg/day 2 . Thus formal uncertainty is only about 4%, and so the estimated 
acceleration is non-zero with very strong statistical significance. This model 
indicates an acceleration starting at 0.024 deg/day 2 at the beginning of the Deep 
Impact data set and reaching 0.028 deg/day 2 by the end of the data set about a 
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day before perihelion. This model was eventually superseded by the acceleration 
profiles described in Sec. 6. 

Hubble Space Telescope and Spitzer Photometry. Photometric 
observations of Tempel 1 taken by the Hubble Space Telescope were obtained 
in 1 997 (Lamy et al. 2001 , 2007) and 2009 (Meech et al. 201 1 ) Each of these 
data sets proved crucial in establishing the spin state in Regions A, B, and C 
because of the dense coverage and high SNR that unambiguously showed the 
light curve morphology. Spitzer flux measurements (Lisse et al. 2005) were 
converted to magnitudes and compared with predictions for visible magnitudes. 
While we chose not to actually fit these data due to uncertainties about the 
relative shape and phasing of light curve extrema, we found that the placement 
of observed Spitzer extrema did agree well with that predicted for visible 
extrema. 

5.2 The Tucson Method. The sidereal spin rates, Sj, are determined separately 
from the rotational phases, Woj. They are determined without reference to model 
lightcurves and, to remove drift in the rotational phase of the observed light curve 
due to changes in the solar phase angle (Surdej and Surdej, 1978; Harris et al., 
1984) and orbital motion, the observing times are adjusted to create a data set 
for each region as it would have been observed by a fictitious observer viewing 
the comet from the direction of the ascending node of the comet at zero solar 
phase angle. After correcting for light time, this ensures that periodicities in the 
modified light curve reflect the sidereal spin of the nucleus. Because of changing 
sub-solar and sub-observer latitude at the elongated nucleus (Fig. 6) the shape 
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of the observed light curve will not be consistent in amplitude, but, because of the 
low obliquity of the comet (1 1 .9°), this effect is expected to be small. We do 
expect, however, that it will introduce some amplitude dispersion in the 
rotationally phased light curves. 

Once the Sj are determined, the rotational phase in each region is then 
estimated by comparing model light curve predictions, based on the Thomas 
shape model and a Hapke phase function, to high S/N Hubble Space Telescope 
data. In the case of Region A, where the Lamy et at. (2001 ) HST data do not 
cover a full rotation period, ground based data are also used in the determination 
of rotational phase. 

Removal of synodic, solar phase angle, and light time effects. The 
observing time for a particular observation is first adjusted to zero solar phase 
angle using a variation of the Phase Angle Bisector (PAB) method of Harris et at. 
(1984). These authors note that as the solar phase angle increases from zero the 
rotational phase of the light curve drifts in the same direction but at approximately 
half the rate. Since the timescale for change in the solar phase angle, a, is much 
greater than the rotational period, the rotational phase Wj(a) of the observed light 
curve at solar phase angle a will be shifted by AWj(a) = p(a/2) from its value at 
zero phase where p is the longitude interval between the meridian that passes 
through the sub-solar point and the meridian that includes the PAB. AW(a) 
corresponds to a timing correction of At(a) « P*AW(a)/360 where P is the 
estimated rotation period. Using the same logic the timing correction for the 
orbital position of the comet at time, t, is At(y) « P*AW s (y)/360 where y is the 
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longitude interval between the meridian that contains the sub-solar point and the 
meridian that contains the ascending node. The total correction to the observing 
time is therefore At = At(a) + At(y) + At(d) where At(d) is the light time correction 
for an observer at a distance d from the nucleus. At varies considerably 
throughout the various observation sets and falls in the range of 0.81 > At > 0.20 
d. 

Preparation of the data for high precision frequency analysis. The data are 
first linearly detrended (this is a small correction as can be seen from Fig.1 1 ) and 
a preliminary value of Sj is determined with the ANOVA (ANalysis Of Variations; 
Schwarzenberg-Czerny, 1996) period-finding algorithm in the commercially 
available Peranso software package (this is available at www.CBABelqium.com ). 
We experimented with the thirteen period finding routines, which include all of the 
major astronomical period-finding algorithms, in the Peranso package on various 
subsets of the data and found that the ANOVA algorithm gave the clearest and 
most consistent results. Other often-used methods, such as FALC (Harris et ai, 
1989) and PDM (Stellingwerf, 1978), gave effectively identical results. 

We have used the uncertainty estimates as calculated in the Peranso 
package. The method used is that described by Schwarzenberg-Czerny (1991) in 
which the uncertainty in the period is taken as the width of the associated peak at 
the mean noise level down from the peak. 

Once the preliminary period estimates have been calculated they are used 
to rotationally phase the data so that they can be examined to estimate bias 
corrections for particular observing runs and remove obviously discordant 
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outliers. In this way the data sets can be “cleaned-up.” This process, removes 
small systematic photometric errors that may arise between groups of 
observations taken under different viewing conditions, with different photometric 
equipment, and under observing conditions of different clarity and seeing. 

The cleaned-up data are then analyzed with the ANOVA algorithm a 
second time to obtain the final estimates of periodicities and uncertainties. This 
process yields our ‘best’ estimate of the period and its uncertainty in each 
observational region. While this two step process is not one that we would 
recommend for finding periodicities ab initio in an arbitrary data set, we have 
confidence in using it here because of the special knowledge produced by the 
Deep Impact encounter photometry. I.e., we know the approximate spin period 
and have a fair idea of the shape of the light curve, knowledge that we use to 
discriminate against spurious period estimates. ANOVA periodograms and 
cleaned-up rotationally phased light curves and for the region A, B, C and the Dl 
data are shown in Figs. 15, and 16 respectively. 

Determination of rotational phases. To determine the rotational phase, 

W 0 j, of the data in each observation region we have relied primarily on HST data, 
which, in regions B and C, provide consistent and well-sampled coverage over a 
complete rotational cycle. The phase is tied into the Thomas et al. (2007) nucleus 
coordinated system by fitting the observed light curve to a model light curve 
computed from the Thomas et al. shape model using the online tool noted earlier. 
The final fits to the data are shown in Fig. 17. The accuracy of fit is estimated at ± 
2° in rotational phase by visual inspection. A detailed examination of Fig. 17 
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shows that there is a distinct possibility of a 180° ambiguity in rotational phase. 
This shown more clearly in Fig. 18. We have used both numerical cross- 
correlation and visual inspection to decide the appropriate alignment of the model 
to the data. However, the results of both methods are, in our opinion, marginal. 
While both of these techniques support the choice of alignment that we advocate 
in this paper the possibility of a -180° error in rotational phase must be 
entertained. 

The fits in Figs. 17 show that the model is delinquent if the range of 
rotational phase associated with the primary minimum of the observed light curve 
and that there is a possible ambiguity of 180° in rotational phase. The final 
estimates for Sj and W 0 j are collected in Table 2 and the evolution of Sj with time 
is shown in Fig. 18. 

Direct determination of the acceleration of the spin rate. Because the comet was 
active during the collection of the Deep Impact data we recognized that it could 
contain direct information on the acceleration of 9P/Tempel 1’s spin. We 
therefore applied the dynamical period estimation methods of Drahus and 
Waniak (2006), which simultaneously yield both a spin rate and its acceleration 
(assumed linear in this case) at the mid-time of the observations (Appendix B 
contains a detailed description of its application to the Deep Impact data). This 
yielded a period of 1 .69961 ± 0.00023 d (21 1 .814 ± 0.029 7d) applicable at JD 
2453530.60510 and a rate of change in the spin rate of +0.020 ± 0.003 °/d/d. As 
can be seen in Fig. 18 this value is consistent with the period change through 
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perihelion that we have found in the Earth based data if the timescale over which 
the torque acts is ~63 days. 

The JPL and Tucson approaches described here are substantially 
independent, and indeed complementary. Through the use of the light curve 
generation tool, the viewing geometry, phase angle correction and light time 
delay are automatically incorporated into the JPL estimate without the careful 
bookkeeping and modification of the time tags required for the Tucson method. 
However, because of its reliance on the light curve generation tool, and in 
particular the Thomas et al. shape model, which did not accurately model the 
light curve at some rotation angles, the JPL approach suffered from somewhat 
poor fits. This led to some irresolvable ambiguities in determining which 
frequency was correct. In contrast, the Tucson approach did not rely on light 
curve models, or any other model, to identify the best fitting rotation rate. In this 
sense the Tucson approach is more robust in determining the rotation rate, while 
the JPL method seamlessly revealed the rotation phase of the comet. 

6. The dynamical evolution of the spin rate of 9P/Tempel 1 

Tables 2 and 3 contain our best estimates of the overall spin state of 9P/Tempel 

1 and the changes that occurred during the perihelion passages in 2000 and 
2005. Table 3 focuses on spin rates and orientation of the polar axis while Table 

2 gives information on rotational phase. The spin rates and acceleration in Table 

3 are the average of the values found in the JPL and Tucson studies. 
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In Figure 18 we plot the observed spin rates and acceleration as a function of 
time. The spin rates derived from Earth-based measurements can be seen to 
systematically increase as the comet passes through succeeding perihelia in 
agreement with the sense, but not the slope, of acceleration that is derived from 
the Deep Impact photometry alone. Referring to the discussion in section 4 we 
assume that the changes in spin rate are the result of a smooth evolution of non- 
gravitational torques through perihelion and from perihelion to perihelion. The 
fact that the spin rate measured by Deep Impact is greater than that measured 
both before and after perihelion and that it is closer in value to that measured 
after perihelion suggests that non-gravitational torques are not symmetric about 
perihelion and that, for a part of the time, positive torques dominate while at other 
times the reverse is true. It seems clear that the dominant effects of non- 
gravitational torques occur well before perihelion passage. 

The acceleration measured during Deep Impact approach allows us to 
compute the magnitude of the torques acting on the nucleus at that time and 
estimate the moment arm that was involved. Assuming a homogeneous mass 
distribution in the interior, a spherical approximation for the shape, and using the 
mass range determined by Richardson et al. (2007), we find that the moment of 
inertia of the nucleus lies between 0.8 - 4.6 x 10 19 kg.m 2 . To achieve the 
observed spin rate acceleration implies that a torque of 0.3 - 2.5 x 10 7 kg.m 2 . s' 2 
was acting on the nucleus. Jewitt (1997) has related the average torque to total 
mass loss in teRMS of a “dimensionless moment arm,” kj, for which he 
conjectured a value of -0.05. With the measurements reported here we can now 
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make an observational estimate of k T for 9P. Schleicher (2007) finds the 
production rate of OH to be ~ 7 x 10 27 mol/s near the relevant time implying a 
water loss rate of only ~2 x 10 2 kg/s. Following Jewitt (1997) we assume a 
characteristic outflow velocity of 10 3 m/s and find 0.005 < k T < 0.04. This 
observational estimate is somewhat lower than Jewitt’s conjecture but the result 
generally substantiates his approach. 

The magnitude of the torque also permits us to say something about the 
stability of the pole direction. If we assume that there is a component of the 
torque acting at right angles to the spin axis for ~60 days (see section 5 above) 
we can estimate the angular displacement under forced precession during this 
time. The angular velocity of precession calculated this way is ~ 1 x 1 0' 9 rad. s' 1 
and leads to a displacement angle of ~ 0.3° during perihelion passage. Unless 
the observed torque turns out to be a very small component of the total torque 
that is operating, which we think is unlikely, we expect that the pole direction of 
9P/Tempel 1 should not precess more than ~1°/perihelion passage. 

Modeling the evolution of the spin state. One of the objectives of this work is to 
predict the rotational state of the nucleus near the time of the Stardust-NExT 
encounter on Feb. 14, 201 1 . To do this a model of the non-gravitational torques 
is required. We have taken two approaches to develop such models - a physical 
approach, in which we attempt to emulate what we know the comets sublimation 
rates and jet structures, and an analytic approach that models the time behavior 
of torques with a prescribed functional form. As with the rotational analysis this 
modeling was done independently at both JPL and in Tucson with similar results. 
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Our physical modeling of torques was based on a heuristic ‘rotating jet’ 
model that has been used earlier to analyze the effects of non-gravitational 
forces on cometary orbits (Chesley and Yeomans, 2005) and in which the 
dependence of torque on heliocentric distance is prescribed by the function g(r) 
(Marsden et ai, 1973). This sublimation “law” has had wide use in non- 
gravitational force studies (Yeomans et ai, 2004a). With this kind of model it is 
possible to match the spin rates and rotational phases of the nucleus across both 
the 2000 and 2005 perihelion passages; however, we also found that it was not 
possible to match the acceleration of the spin rate measured on Deep Impact 
approach at the same time. For this reason modeling based on jet torques was 
abandoned by both groups. An alternative physical approach is to consider 
torques associated with widespread sublimation of H 2 O over the surface of the 
nucleus employing the Thomas et ai (2007) shape model. Such a study is 
already underway and Samarasinha (2010, private communication) reports that 
initial results indicate that model period changes through perihelion are 
somewhat similar to those measured in this work and that both phases of positive 
and negative net torque can occur during perihelion passage. Since this work 
was already in an advanced state we decided not to duplicate the effort and 
focused instead on analytic modeling. Since the analytic models produced at JPL 
and Tucson are quite different we present them separately. Even though the 
methods used by the two groups are quite different we shall find that they yield 
very similar descriptions of how the net torques currently operate during 
perihelion passage: on approach the spin rate first decreases, passes through a 
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minimum, and then accelerates rapidly through perihelion. After perihelion the 
spin rate goes through a maximum and then decreases to a stable level as the 
nucleus moves away from the sun. 

The JPL torque model. Given the solutions detailed in Table 2 we can try to link 
together the various quiescent periods by examining the rotation phase runoffs 
induced by the torques encountered during the active intervals. We start with the 
2005 perihelion passage since this case affords information on the comet’s 
rotation state near the time of perihelion, as well as before and after. From fits to 
the Deep Impact approach data we know that the rotation phase W=219° on 
2005 Jul 5.0 at perihelion. Meanwhile the pre- and post-perihelion quiescent 
solutions (Regions B and C in Table 2) predict W= 225° and W=104°, 
respectively, at perihelion. These imply that the accelerated comet gained 354° 
or -6° in rotational phase as it approached relative to a hypothetical 
unaccelerated comet. We call this the runoff. Similarly, but working in reverse 
time, the post-perihelion runoff is 245° relative to the unaccelerated post- 
perihelion solution. Thus according to these estimates the combined pre- and 
post-perihelion runoff is 239°(modulo 360°). 

Using a similar approach at the 2000 perihelion, but without any active 
period constraints, we find that the combined runoff should have been about 
208°, again modulo 360°. 

We could not identify a torque profile that meets all of the constraints and 
leads to 354° runoff at perihelion, while we have developed an ad hoc model that 
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does lead to -6° runoff (Fig. 19). This model assumes piece-wise constant 
accelerations, and thus piece-wise linear spin rates. The key feature is that there 
is a modest (and so far not directly observed) deceleration before the dramatic 
acceleration seen in the Deep Impact data begins. This deceleration period 
would be about -0.0013 deg/day 2 and would start 400 days pre-perihelion and 
end 660 days post-perihelion, which corresponds well to the active periods seen 
in Fig. 4. Across this background deceleration is superimposed an acceleration of 
approximately 0.029 deg/day 2 that acts only for the 86 days before perihelion. 
This scheme would arise from a diffuse negative torque that acts for most of the 
active period and a single, strong, seasonal, jet-like active region that only acts 
for the three months prior to perihelion. 

While this model is obviously crude and lacks a detailed physical basis, it 
does meet the observational constraints and serves one of the key purposes of 
the project, which is to predict the spin state of the comet at the epoch of the 
Stardust-NExT flyby in mid-February 2011. Here we can assume that the comet 
will essentially do the same as it did in 2005, or we can assume that there is 
some secular change and extrapolate from 2000 to 2005 to 201 1 . The former 
assumes implicitly that the perihelion-to-perihelion changes are best modeled as 
a random walk, while the latter would be most appropriate under the assumption 
that the comet is changing in a predictable way. Given the paucity of information 
we have about the variability of comet nucleus activity, each of these 
perspectives are equally defensible. For the present model we assume the 
random walk hypothesis and thus predict that the 201 1 runoff will be -6° to 
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perihelion (201 1 Jan 12.4). The additional 33 days from perihelion to the 
Stardust-NExT encounter (about 201 1 Feb 15) should accumulate an additional 
55° of runoff, amounting to a total of roughly 49°. The quiescent, post-2005 
solution predicts W = 349° at 201 1 Feb 1 5. Adding in the presumed runoff of 49° 
yields the accelerated prediction, W = 38°. The spin rate at that time should be 
213.53°/day giving W 0 = 172° and, on 201 1 Feb 15.0, the sub-solar longitude is 
calculated to be 328°. The mission target is W = 98° at closest approach and so 
we should delay the flyby by 60°, i.e. delay the encounter to 201 1 Feb 1 5.28. 
This is a 10 hour delay from the nominal (i.e., before February 2010) arrival time 
of 2011 Feb. 14.8625. 

The Tucson or “Gauss” analytic model. Here we divide the component net 
torque, T(t), at time t in Eq. 4 into two parts, to give the model the flexibility to 
represent both negative torques and positive torques separately as follows: 

T(t) = T 1 *exp((t-t max1 )/x 1 ) 2 ) + T 2 *exp((t-t max2 )/x 2 ) 2 ) (9) 

The Gaussian shape assumed here is quite arbitrary and was chosen simply 
because it had the right general character for torques generated by sublimation 
and, more importantly, it allowed the integration for the spin rate to converge. 
The shape should be thought of as an interpolation function. We did experiment 
with Lorentzian shapes, but found that the extended wings gave unrealistic 
torques that extended too far into the region around aphelion. 
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The model therefore has six parameters, Ti,T 2 , t max i, t max2 , t-i, and x 2 and, 
when applied to the 2005 perihelion passage, seven constraints - the observed 
spin rate and rotational phase at three times and the acceleration of the spin rate 
at a single time. A MathCad program was developed to solve Eq. 4 as a function 
of assumed values for the parameters t max i, xi, and x 2 . Preliminary values for T i , 
T 2 and t max2 were first estimated using the three observed spin rates and the 
acceleration at Deep Impact encounter. Given the initial rotational phase in a 
torque free region, Eq. 6 could then be integrated to provide an estimate of the 
rotational phase spin rate at the time of the other two observations. The 
parameters t max i, xi, and x 2 were then adjusted (by hand) in an iterative cycle until 
the model rotational phases agreed with the observed values to better than 1°. 

At the 2000 perihelion we have only four constraints and it is not possible 
to solve for all of the parameters. We have therefore assumed that the 
parameters that depend primarily on the geometry of the perihelion passage, 
tmaxi, t max2 , xi, and x 2 , have the same values as in 2005. We expect this to be a 
reasonable assumption providing that any precession of the spin pole (or the 
orbit) is, as we have assumed, negligible. In this way it is possible to calculate 
appropriate values of Ti and T 2 for the 2000 perihelion passage. 

The results of these calculations are collected in Table 4 and the torque 
profile and spin rate evolution for the 2005 perihelion passage are shown in 
Figure 20. We see that the observations, as expected, require T-i and T 2 to have 
opposite signs. The predicted net torque profile initially moves to a negative 
value before achieving a positive maximum some 32 days before perihelion. In 
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addition, as expected from the decrease in the amount of period change through 
the 2005 perihelion relative to that at the 2000 perihelion, the magnitudes of the 
model torques are also seen to decrease from 2000 to 2005. 

We use the parameters determined at the 2000 and 2005 perihelion 
passages to predict what the torque profiles may look like in the 201 1 perihelion 
passage. First, we again assume that the parameters that depend primarily on 
geometry remain the same in 201 1 as in 2005. Second we linearly extrapolate 
the values of T i and T 2 found in 2000 and 2005 to 201 1 . The resulting values of 
the Ti and T 2 are shown in Table 4 and the predicted spin rate profile for 201 1 in 
Figure 21 . Putting these three solutions together we can plot our estimate for the 
complex evolution of the spin rate throughout the period 1997 - 2010 and the 
prediction for 201 1 . This is done in Figure 22. 

7. Discussion and predictions for Stardust-NExT encounter 

The observations collected during the Deep Impact World-wide campaign cover 
an interval of thirteen years and two perihelion passages and clearly imply (Fig. 
18) a roughly “stepwise” increase of the spin rate of 9P/Tempel 1 as the comet 
passes through succeeding perihelia. In addition, most of the torque must have 
been applied well before perihelion. The spin rate and its acceleration measured 
from the Deep Impact approach photometry imply even greater complexity in 
requiring that phases of both positive and negative net torque operate through 
the observed perihelion passages. The observations and modeling appear to 
support the idea that the torque producing outflow is dominated by the 
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sublimation of H 2 O over a large fraction of the surface of the nucleus and is a 
result of its shape. Modeling indicates that the production of net torque is unlikely 
to be dominated by observed jet structures except, possibly, in the period just 
before perihelion passage. While the spin pole direction is certain to precess as a 
result of these torques, the rate is expected to be small, < 1°/perihelion passage. 
If it becomes possible to measure the displacement of the direction of the spin 
pole at the Stardust-NExT encounter in 201 1 we should be able to deduce more 
about the distribution of the outflow with respect the surface of the Thomas et al. 
(2007) shape model. 

Relationship to the water production rate. Observations of OH production by Osip 
et al. (1992), Schleicher (2007), and Cochrane et al (2009) and, in particular, 
their interpretation by Schleicher, indicate that the water production rate peaks 30 
- 60 d before perihelion and has moved closer to perihelion by 5-10 d between 
the 1983 and 2005 perihelion passages. Apportioning this change equally 
among the intervening perihelia this gives roughly a shift of 1 - 3 d/perihelion. In 
addition, the production rate decreased by 42% between 1983 and 2005. 

These changes are qualitatively consistent with the present observations in that 
the period change, and therefore the implied production rate of the operating 
outflow, decreased between 2000 and 2005. Also the (model) torque peaks well 
before perihelion (32 d) as does the observed H 2 0 production rate. The H 2 0 
production rate implied by the model at the peak torque fell by 5% between 2000 
and 2005. 
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Our observations and analytic modeling do not provide a quantitative 
explanation for the shift of the H 2 O production peak by 1 - 3 d/perihelion 
passage. But we presume that this is the result of small changes in the 
distribution of sublimation over the nucleus surface from perihelion to perihelion 
or, possibly, slow precession of the spin axis. 

Predictions for the Stardust-NExT encounter with 9P/Tempel 1. The 
primary objective of the mission is to see what changes occurred on the regions 
of the surface previously imaged in the Deep Impact mission as a result of the 
activity during a single perihelion passage and, secondarily, to characterize the 
artificial crater formed by the impactor spacecraft. In January, 2010, when time- 
of-arrival trajectory maneuvers were being planned, the spacecraft was projected 
to arrive at the comet on UT 2011, February 14 at 20:42:34.8 (JD 
2455607.36290) some 34 d after perihelion passage. To have the best chances 
of accomplishing these objectives the spacecraft should arrive at encounter when 
at least 25% of these regions and the crater are in daylight and visible to the 
spacecraft camera system. Figure 23 is a contour plot which shows the 
encounter conditions that must prevail in order to achieve the above objectives. 
By providing a time-of-arrival (which sets the sub-solar longitude at encounter 
time) and a spacecraft B-plane angle that avoids the dark blue area it is possible 
to ensure that the mission will achieve its objective. 

In the JPL study, in order to make a prediction for the spin state at 
Stardust-NExT encounter, the spin evolution model that was developed for the 
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2005 perihelion passage was simply transferred to the 201 1 time frame by 
modifying the initial spin rate and rotational phase to that appropriate for region 
C. This approach was based on the assumption that changes in the torque profile 
from perihelion to perihelion cannot be easily predicted and that the comet is 
likely to behave at one perihelion essentially as it did at the previous perihelion. 
This yielded a spin rate of 213.52 °/d and a sub-solar longitude of 328°W at 
encounter on UT 201 1 Feb 15.0 (299°W on Feb 14 20:40). 

In the Tucson study, the value of the parameters Ti and T 2 in the Gauss 
model that are appropriate for the 201 1 perihelion passage are obtained by a 
linear extrapolation of their values in 2000 and 2005. The values of the other 
parameters in the Gauss model are, as explained earlier, taken as the same as 
in 2005. The model yields a spin rate of 213.47 °/d, W 0 = 15.9°, and a sub-solar 
longitude of 242 W on 201 1 Feb 14 20:40. The two models differ by 57° in their 
prediction of sub-solar longitude. In order to be specific, and since the two results 
are reasonably close (i.e., both well within a single quadrant), we simply take 
their average as our final estimate for the sub-solar longitude at nominal 
encounter and treat the difference an indicator of the level of uncertainty in the 
result, i.e., we take the sub-solar longitude at the nominal time of encounter (UT 
201 1 Feb 14 20:40) to be 271 ± 29° W. This result is plotted on Figure 23 as a 
black dot and bar and shows that a trajectory correction maneuver to delay the 
time of arrival by at least 8h is needed to be sure of attaining the primary science 
objective. In Figure 24 we show, using the Thomas et al. (2007) shape model, 
the predicted aspect of the nucleus as seen from the Stardust-NExT spacecraft 
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at the nominal encounter time (LHS) and the case where the encounter time has 
been delayed by 8h (RHS) where there is not only ample Deep Impact terrain to 
observe but also what should be a spectacular view of the artificial crater. A burn 
to accomplish a ~8h delay in arrival time (the maximum allowed by the available 
fuel) was performed in February, 2010. 

Finally, we note that there is an important caveat to the above prediction: 
there is, as noted in section 5, the possibility of an error of 180° in the predicted 
rotational phase. If this is the case then the nominal encounter time will achieve 
excellent viewing of Deep Impact terrain at a sub-solar longitude of -91° W 
without a trajectory correction maneuver. With a time delay of 8 h, implementing 
the trajectory correction moves the sub-solar longitude at encounter to ~162±25° 
W. As can be seen from Figure 23 this still allows the possibility that the primary 
science objective will be achieved, but that imaging of the artificial crater would 
be unlikely. 

8. Conclusions. 

In this paper we have provided a detailed analysis of light curve information from 
the Deep Impact World-wide observing campaign and data obtained from the 
Hubble and Spitzer Space Telescopes and the Deep Impact mission that were 
obtained between 1997 and 2009. This analysis shows: 

1 . The spin rates (periods) changed in an approximately stepwise manner 
through the 2000 and 2005 perihelion passages. From 209.023±0.025 °/d 
(1.7223 ± 0.0002 d) prior to 2000 to 210.448 ±0.016 °/d (1.7106 ± 0.0001 


47 



d) between 2000 and 2005 and then to 21 1 ,814±0.029 7d (1 .6996 ± 
0.0002 d) during the Deep Impact approach to 21 1 ,625±0.012 7d (1 .701 1 
± 0.0001 d) post-2005. 

2. The period shortened by 16.8 ± 0.3 min during the 2000 perihelion 
passage and by 13.8 ± 0.2 min during the 2005 perihelion passage. 

3. The angular acceleration was 0.024 ± 0.003 °/d 2 during the Deep Impact 
approach. 

4. In 2005 the angular acceleration was not symmetric about perihelion and 
most occurred well before perihelion passage. 

5. Sublimation outflow of H 2 O over most of the surface of the elongated 
nucleus with a possible contribution from jet activity just prior to perihelion 
are likely the causes of the torque. 

6. The level of torque required to explain the Deep Impact observations 
suggests that precession of the spin axis is small, i.e., < 1°/perihelion 
passage. 

7. The trend in the net change of spin rate through the two perihelion 
passages is in a direction expected from the published trend in H 2 O 
production rates observed at the 1983 through 2005 perihelia. 

8. The observed peak in H 2 0 production rate some 30 - 60 days before 
perihelion is in concert with the predicted peak in torque at 32d before 
perihelion. The H 2 0 production rate implied by the Tucson torque model 
fell by 5% between the 2000 and 2005 perihelia. 
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9 . Analytic models have been used to predict the rotational state at the 
comet at the planned Stardust-NExT mission encounter in 201 1 . On UT 
Feb 14, 2010 20:40, the nominal encounter time in January, 2010, we 
predict that the sub-solar longitude on the nucleus will be 271 ± 29° West 
longitude (the average of the JPL (299°W) and Tucson (242°W) 
determinations). The spin rate is predicted to be 213.5 ± 0.2 °/d. 
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Table 1: 


The data used in the rotational analysis 


Region 

Number of 
observations 

Time span 
(days) 

Standard 

deviation 1 

Heliocentric 2 
range (AU) 

Dates 

A 

307 

600 

±11 

4.5 -2.2 

11/23/97-7/15/99 

B 

541 

1000 

±12 

4.1 -3.5 

8/16/01-5/11/04 

Dl 

2888 

50 

±23 

1.6- 1.5 

5/15/05-7/3/05 

C 

1074 

1194 

±13 

3.6 -3.3 

9/24/06-12/30/09 


1 In brightness units; the mean nucleus brightness = 100 units. The large 
standard deviation in the Deep Impact data is due to noise from the coma that 
was present in the signal. 2 For regions A, B, C the heliocentric range passes 
through aphelion at 4.74 AU. 
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Table 2: 


Results of rotational analysis: Values of W(t)j = Woj + Sj*At for each 
region. Units: W 0 j are in degrees, Sj are in °/d, and dS D i/dt is in deg/d/d. 
At = (t(JD) - 2451545.0); ^he JPL study assumed the Tucson value 
for S A was correct. The spin rates for the Deep Impact data are given 
for the case of zero acceleration. When evaluating the values of Woj 
between the two studies note that part of the difference is due errors in 
Sj propagating back to the reference time. The intrinsic lightcurve fitting 
error in Woj (~ ± 2°) is therefore much less than the differences would 
imply. Formal errors are given in the text. 


JPL method 

Tucson method 

W(t) A = 69 + 209.023 ± 0.025 u At 

W(t) A = 72 + 209.023 ± 0.025*At 

W(t) B = 289 + 210.438 ± 0.0xx*At 

W(t) B = 280 + 210.458 ± 0.01 6*At 

W(t) D i = 255 + 21 1 .849 ± 0.0xx*At 

W(t) D i = 299 + 21 1 .862 ± 0.030*At 

W(t) c = 299 + 21 1 .626 ± 0.0xx*At 

W(t) c = 301 + 21 1.623 ± 0.01 0*At 

dSoi/dt = 0.024 - 0.028 

dS D |/dt = 0.020 ± 0.003 
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Table 3. 


The observed spin state of 9P/Tempel 1 and its changes through two perihelion passages 
(CY 2000 and 2005). * This estimate includes the effect of acceleration (assumed linear) 
during Deep Impact (Dl) approach. + This estimate assumes zero acceleration in the spin 
rate during Deep Impact ( Dl) approach. Information on rotational phase is included in 
Table 2. 


Orientation and sense of rotation: 

References 

Sense of spin 

Direct 

Thomas et at. (2007) 

Direction of pole (Positive or 
North) 

RA=294°, Dec=73° 

(± 5° on the sky; J2000) 

Thomas et at. (2007) 

Motion of the pole 

Precession possible? Change not detected over 8 apparitions. 

< 1°/perihelion passage. Nucleus in essentially fully relaxed SAM 
(Short Axis Mode) state of rotation. 

Yeomans et al. (2004b); 
Schleicher (2007); This 
paper 

1 1 

Rotation rates: 

Epoch (days) 
(JD -2450000) 

Period 

(days) 

Angular rate (°/day) 


Pre - 2000 

1224.88 

1.7223 ±0.0002 

209.023 ± 0.025 

This paper 

2001 - 2004 

2477.37 

1.7106 ±0.0001 

210.448 ±0.016 

This paper 

Deep Impact Approach* 

3530.60 

1.6996 ±0.0002 

211.814 ±0.029 

This paper 

Deep Impact Approach* 

3544.50 

1.6993 ±0.0002 

21 1.856 ±0.030 

This paper 

Post - 2005 

4297.00 

1.7011 ±0.0001 

211.625 ±0.012 

This paper 

1 1 

Rotational acceleration: 


dP/dt (days/day) 

doj/dt (deg/da/) 


Deep Impact approach 

3530.60 

-1.9 ±0.1 x 10' 4 

0.024±0.001 

This paper 

1 1 

Net change in rotation: 


AP(minutes) 

Aw(7day) 


2000 perihelion passage 

1546.13 

-16.8 ±0.3 

1.425 ±0.030 

This paper 

2005 perihelion passage 

3556.82 

-13.7 ±0.2 

1.177 ±0.020 

This paper 
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Table 4. 

Parameters for the Tucson “Gauss” torque model. The parameters are 
defined in Eq. 9. Details of the fitting process to the 2005 and 2000 
perihelion spin rate and acceleration data is described in the text. The 
units of Ti and T 2 are °/d/d. 


Perihelion 

Ti 

t 2 

Ti (d) 

^ 2 (d) 

tmaxl(d) 

tmax2(d) 

2000 

-0.112612 

0.132565 

180 

159 

-16.9 

-22 

2005 

-0.110612 

0.129355 

180 

159 

-16.9 

-22 

2011 

- 0.108612 

0.126145 

180 

159 

-16.9 

-22 
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Appendix A: 

Description of the Deep Impact MRI Photometric Reductions 
(Contributed by Fabienne A. Bastien) 


Photometric measurements of comet 9P/Tempel 1 were performed on images 
taken with the Medium Resolution Instrument on the Deep Impact spacecraft 
during the approach phase from 1 May 2005 to approximately 6 hours before 
impact on 4 July 2005. The MRI is a 12 cm aperture Cassegrain telescope with 
a 2.1 m focal length. Images were taken through clear filters that have a center 
wavelength of 650 nm and are uncoated and not band limited (Hampton et al., 
2005). 

A total of 3014 images were analyzed: 595 science images and 2419 
optical navigation images. The measurements were based on circular apertures 
ranging from 5 to 30 pixels in diameter with nucleus at the center of the aperture. 
The data consist of science images taken with the clear 1 and clear 6 filters and 
optical navigation images taken with the clear 1 filter. In the following discussion 
the images are displayed with lines increasing up and with samples to the right. 
Figure Al (taken from Klaasen et al., 2008) shows a full frame image with the 
quadrant nomenclature used below. 

Properties of the Science data : Photometric measurements were made from the 
reversibly calibrated ("RADREV") science images. These images have had the 
standard pipeline corrections applied to them: bias and dark frame subtraction, 
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flat-field corrections, etc. They have not, however, been "cleaned" to remove 
artifacts such as cosmic rays. All images were taken in one of two sub-frame 
modes: 256x256 pixels for most of the approach sequence, and 512x512 pixels 
for the last 1 .7 days of approach. The images display a number of problems not 
accounted for in the standard reduction. All images were affected by a horizontal 
striping of a few DN in amplitude caused by electrical interference. Additionally, 
the bias levels for these images are only known to the nearest full DN, and the 
bias level for each quadrant is different. Not correcting the images for this 
interference and, more importantly, for the imprecise bias subtraction can 
introduce systematic errors in the photometry in which the early approach images 
are primarily affected. 

A further problem concerns the pixels in the two rows surrounding the 
horizontal boundary between the upper and lower halves of the CCD that are 
each 1/6 of a pixel smaller than the other pixels of the CCD due to the way the 
readout clocking was designed. This increases the point spread function of 
objects that overlap this boundary by 1/3 of a pixel. Flux measurements 
therefore tend to be greater at the boundary because the flat fielding during the 
pipeline processing assumes that all pixels are the same size. In approximately 
80% of the images, the centroid of the nucleus lies within 20 pixels of this region; 
images taken through the clear 6 filter are particularly affected. Thus, most of the 
photometric measurements need to be corrected for this effect. 

Properties of the navigation data. The raw optical navigation images are in a 
different format from the science images. Each image consists of a number of 
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square "snippets", each less than 400x400 pixels in size. Each "snippet" is 
centered on an object deemed interesting by the navigation team, usually a star 
or 9P/Tempel 1 . Each set of snippets is integrated into a single image that is 
1008x1008 pixels in size; the navigation data have neither serial-overclock nor 
parallel-overclock pixels. Therefore, quadrant bias must be determined via other 
methods. For a more detailed description of the navigation data, please see the 
Deep Impact Navigation Images Report included in the Deep Impact 
documentation data set, DI-C-HRII/HRIV/MRI/ITS-6-DOC-SET-V1 .0, which is 
archived by the Planetary Data System. 

Optical navigation data suffers from the same striping noise as the science 
data. However, because these data are in a different format from the science 
data, a different algorithm had to be applied to correct it. This procedure was 
also used to remove the bias. For the navigation data, the centroid position of the 
comet was sufficiently far away from the horizontal boundary between the upper 
and lower halves of the CCD (more than 30 pixels, on average) so that the 
photometric measurements were not affected by the smaller size of the pixels 
there. 

Photometric reduction of the science data. The general procedure for the 
analysis of the science data is as follows: We start with RADREV calibrated data 
and then remove the horizontal striping. We then convert the data back to DN/s 
and perform the aperture photometry. These data are then corrected for the 
quadrant boundary effect (where needed). In order to remove the horizontal 
striping, the image is first divided into two halves: quadrants B and D to one side 
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and quadrants A and C to the other. The process is similar for each half: a region 
is defined, avoiding the overclock pixels, which will be used to determine the 
background. This is a two-dimensional array 25 pixels wide and whose length 
depends on the size of the image (256 pixels or 512 pixels). We take a resistant 
mean across each row of this array and store into a new array that is 1 pixel wide 
and 256 pixels or 51 2 pixels long. This then subtracted from each column of the 
half-image. Once the background noise is removed, circular aperture photometry 
is performed with apertures ranging in size from 5 pixels to 30 pixels in diameter. 
If the aperture falls across the boundary between the upper and lower halves of 
the CCD, the following procedure is used to correct the flux: We create a sub- 
image centered on the comet that is slightly larger than the aperture. This sub- 
image is then divided into two parts: part 1 is the portion of the image located 
above the quadrant boundary and part 2 corresponds to the part below the 
boundary. We measure the flux contained within the original aperture (i.e. with 
the center of the aperture at the original centroid position) that is in part 1 . The 
centroid of the aperture is then shifted up by 1/6 of a pixel, and the flux re- 
measured. The difference between these two measurements is half of the 
necessary correction to be added to the total flux. The procedure is repeated for 
part 2, only this time the centroid of the aperture is shifted down by 1/6 of a pixel. 
Finally, we take the two central rows of the CCD (rows 51 1 and 512), and 
measure the flux from these two rows that contained within the aperture. One 
third of this value is subtracted from the total flux. Figure A3 illustrates this 
procedure. 
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Photometric reduction of the navigation data : The process is similar to that used 
for the science images except that the correction for the quadrant boundary is not 
applied. The background removal algorithm that we apply to the navigation 
images differs from that used for the science images in that the edge of the 
relevant snippet is used estimate the background instead of the edge of the 
image. Since the apertures used do not cross quadrant boundaries, it is only 
necessary to use one edge of the snippet. This procedure simultaneously 
determines the value of the background in the snippet and the bias. 

Unfortunately this procedure does not work for images that were taken 
during the last week of approach. At this point, the comet's coma contaminates 
the entire snippet centered on the comet, and rarely are there any other snippets 
within the same quadrant (this is always quadrant D during this time period). In 
this case the value of the background is measured from the raw science images 
taken closest in time to the navigation image under consideration. A slightly 
different bias value is subtracted from the navigation images (358.5 DN) than 
from the science images (359 DN); however since the bias values applied to the 
science data are only determined to the nearest full DN, these two numbers are 
consistent with one another. Note that horizontal striping is not removed from 
images taken during this time frame. At this point in the sequence the comet is 
bright enough for the effect to be negligible. In 12 cases the comet was close 
enough to the quadrant boundary to require correction to the flux. Since these 
images comprise less than 0.5% of the total number of navigation images, these 
data were simply omitted from the analysis. 
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Appendix B 

Investigation of Deep Impact photometry with the dynamical techniques. 

Contributed by Michal Drahus 


Deep Impact approach photometry (Section 3) was obtained close to perihelion, 
when the comet was active - thus providing an opportunity to seek angular 
acceleration in this data set alone. This would manifest itself as a small deviation 
from the constant periodicity. We performed such an analysis using the 
dynamical techniques introduced by Drahus & Waniak (2006). 

We used their dynamical implementations of two classical algorithms: the 
Phase Dispersion Minimization (PDM), introduced by Stellingwerf (1978) and 
later improved by Drahus & Waniak (2006) to weight input data points according 
to their errors (which they called DPDM), and a Least Squares fit (hereafter 
harmonics fit) of a sum of harmonics, which also weights the input data. Both 
methods return a variance ratio R, which is a function of the rotation frequency 
coo and acceleration doo/df, and whose minima indicate the best dynamical 
solutions. Although the techniques allow for any a-priori law controlling the 
accelerating torque, in our analysis we assumed it was constant. This implies a 
constant doo/d t and a linear evolution of go with time. When the frequency is 
labelled with index zero, it refers to a specific moment of time f 0 ; throughout this 
section to is the middle moment of the Dl data set, which is June 9, 2005, 
2:31:20.6 UT (JD2453530.60510). 
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Discovery of the secular spin-up of Tempel 1 made it possible to predict 
the frequency ooo and acceleration doo/df for the moment to, and consequently to 
limit our analysis to the vicinity of the expected solution. Hence we investigated 
too between 0.024 and 0.025 h' 1 (which corresponds to the periods Po between 
1.667 and 1.736 d), and doo/df between -5.0 and +5.0 * 10' 7 h' 2 . The DPDM was 
used with 20 to 200 bins (and always 5 covers), and the harmonics fit with 3 to 7 
harmonics (including the base sinusoid). However, the results were found to be 
very weakly dependent on the settings. For the sake of clarity we present the 
solutions from the DPDM with 60 bins and from the harmonics fit with 5 
harmonics (Fig. B1 ), and adopt their mean value as the final dynamical solution. 
The solution is oo 0 = 0.0245155 ± 0.0000033 h' 1 (P 0 = 1.69961 ± 0.00023 d) and 
doo/df = +0.97 ± 0.15 x 10' 7 h' 2 , which unambiguously confiRMS slow spin-up of 
the nucleus. The parameters are significantly correlated, with the correlation 
coefficient of -0.70, which is a consequence of non-uniform distribution of the 
input data points. Note, that although data phasing is influenced simultaneously 
by ooo and doo/d t, for uniformly distributed points analysed with respect to the 
middle moment t 0 , the correlation would be weak or completely removed. 

Errors and the correlation coefficient were estimated following the Monte-Carlo 
approach of Drahus & Waniak (2006). We take as the noiseless reference light 
curve the harmonics fit with 5 harmonics, calculated separately for the solutions 
from both methods, we simulated 1000 realizations of noise for each of them, 
and determined the covariance matrices. Consequently, the mean matrix 
provided the errors and the correlation coefficient for the mean solution. It is 
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worth noticing, that the dispersion of both individual solutions about the mean is 
only 15% of the error in frequency and 56% of the error in acceleration, which 
shows an excellent consistency of both algorithms and suggests that 
uncertainties are reliably estimated. 
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Figure captions 


Figure 1. 9P/Tempel 1 R(1 ,1 ,a) magnitudes for 1997 through 2010 from the data 
of Meech et al. (2011) and Lamy et al. (2001 , 2007). The ground based data 
are filled circles and the HST data are in red. The data set contains a subset 
of published V magnitudes that have been converted to R using (V-R) = 0.50. 
The short black vertical lines show the times of opposition and the dashed 
vertical lines denote the time of perihelion passage. The horizontal orange 
bars show the range of data in each region that is used in the rotational 
analysis. Note how the comet brightens as the solar phase angle decreases 
near opposition. 

Figure 2. Phase laws. The R(1 ,1 ,a) data (black dots) with various phase 

functions proposed for 9P/Tempel 1 . Shown are the integrated phase function 
of Li et al (2007) (dotted line), the phase law polynomial of Belton et al (2005) 
(dashed line), and a linear correlation to the data (continuous line). We use 
the latter (R(1 ,1 ,a) = 14.905 + 0.0449a) to correct the magnitudes to zero 
phase angle (a). 

Figure 3. R(1 ,1 ,0) versus calendar year. The magnitude data has been corrected 
to zero phase angle using a linear phase law of 0.0449a mag/deg. The flat 
trend around the time of each aphelion passage indicates that coma there is 
negligible and the scatter in that region is almost entirely due to the spin of 
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the nucleus. The red circles denote HST data. The horizontal orange bars 
shows the range of data in each region that was used in the rotational 
analysis. 

Figure 4. R(1 , 1 ,0) versus time from perihelion passage. The filled points are 
relative to the Jan 2, 2000, perihelion passage and the open circles are 
relative to the July 5, 2005, perihelion passage. The open triangles were 
taken in the fall of 2008 and refer to the 201 1 perihelion. This figure shows 
the repeatability of the data from one perihelion to the next and also the 
asymmetry of the light curve about perihelion (see also Figure 5). 

Figure 5. R(1 ,1 ,0) versus heliocentric distance. The brightness asymmetry about 
perihelion is clearly evident. The filled points are pre-perihelion, the open 
circles are post-perihelion. Notice that the mean trend becomes independent 
of heliocentric distance beyond 4 AU, indicating that light scattered from the 
nucleus dominates the signal. 

Figure 6. Sub-solar (solid line) and sub-earth (dashed line) latitudes (left 
ordinate) for 9P/Tempel 1 from 1997 to 201 1 . The data, R(1,1 ,a), are also 
shown (right ordinate). The sub-earth latitude varies from -5 deg at the 2002 
opposition to -14 deg at the 2004 opposition and this change of viewing 
geometry has an effect on the shape of the light curve. The increase in sub- 
earth and sun latitudes from negative to positive latitudes just prior to 
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perihelion passage is important in the interpretation of H 2 O productions rates 
(Schleicher, 2007) and consequently the torques on the nucleus (see Section 
6 ). 

Figure 7. The orbits of 9P/Tempel 1 and the Earth projected onto the ecliptic 
(since the comet inclination is low, its true orbit rises above the plane of the 
figure by, at most, 0.482 AU). Shown are the directions of the comet from the 
sun at each opposition covered by the data. The wide spread of these 
directions shows why timing corrections must be made to the observations if 
the data from different oppositions are to be phased together and a sidereal 
rotation period determined. Also shown (dashed line) is the projection of the 
direction from the ascending node of the nucleus equatorial system (X = 
202.03, p = 9.24) toward the nucleus that is used as a basic reference 
direction when making corrections for synodic and illumination effects. 

Figure 8. Fraction of signal due to the nucleus in each data point of the Deep 
Impact approach photometry. Evidently coma dominates the photometric 
signal for most of the data until about 4 days before impact (see text for a 
detailed explanation). 

Figure 9. Deep Impact approach photometry. The top two panels show a 
normalized version of the science and navigation data separately after the 
clean-up process that removed variability due to mini-outbursts (Farnham et 
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a/., 2007; Belton et al., 2008). The data, in ‘data numbers’ (DN), are 
normalized by the factor r 2 d 2 /C n f(a) as described in the text. Also shown are 
4 th -order, least-square fit, polynomials fitted to the data. These represent the 
general run of the sum of the underlying coma signal and the mean nucleus 
signal. In the lower panel is the difference between the data and the 
polynomials, which is taken to represent the rotational variability of the 
nucleus. 

Figure 10. Plots of post- and pre- 2005 perihelion brightness in relative 
brightness units (mean nucleus brightness = 100 units) after the mean 
brightness of the nucleus has been removed. The dashed horizontal lines 
show the mean brightness of the nucleus for comparison purposes. These 
plots make it easy to evaluate the relative contributions of the coma and 
nucleus to the total brightness and determine the heliocentric distances 
between which the production of a coma is insignificant, i.e., from 4.1 AU 
post-perihelion to 3.2 AU pre-perihelion. The curves are polynomials, which, 
as can seen, only approximately represent the run of the coma levels. They 
are not used in the analysis and are shown only to indicate the general trend 
of coma brightness. 

Figure 1 1 . Brightness data for 9P/Tempel 1 used in the rotational analysis. The 
data are shown after the clean-up process (see text) and displayed as a 
function of observing time. The three observing regions defined in Figure 1 
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are shown separately. The mean brightness of the nucleus (100 units) has 
been removed. 

Figure 1 2 - An example of the fit of data in the JPL analysis over a limited 
interval around an assumed to (in this case in Region A) to the predicted 
model light curve (black line). 

Figure 13 - An example of an alternate method to determine the spin rate used 
in the JPL analysis where the RMS of residuals are minimized over fits to a 
short, dense, data sample over a range of assumed spin rates (co). 

Figure 14 - Model light curve fits to the Deep Impact approach photometry in the 
JPL analysis. The dark black line is for the Hapke phase function. The colors 
of the data points represent different regions of the fit that were used to 
estimate the acceleration of the spin rate. 

Figure 1 5 - ANOVA periodograms for the data in Figs. 9 and 1 1 . 

Figure 16 - The cleaned-up data for each observation region phased with the 
spin rates in Table 3. 

Figure 17 - Fits of model light curves to the HST data that determine the 
rotational phase in each observation region. In the top two panels we show 
the preferred fit based on numerical cross-correlation and visual inspection. 
Both of these methods give only marginal assurance of the best fit (i.e., a 
180° ambiguity is possible). Note the discrepancy in the fit near the primary 
minimum in both sets of data. In the bottom two panels we show the fit of the 
model to a combination of Region A HST and ground-based data again 
illustrating the possibility of a 180° ambiguity. 
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Figure 18. Spin rate results for 9P/Tempel 1 . The black points are the spin rates 
deduced from the data in regions A, B, Deep Impact, and C and listed in 
Table 3. The dashed line indicates the slope determined from the Deep 
Impact approach photometry. The horizontal red lines show the time period 
covered by the data in each region that was used to estimate the spin rate. 
The predicted spin rate (212.522 ±0.174 °/d) near aphelion following the 
201 1 Stardust-NExT encounter is shown in aqua. The observational error 
associated with the observed points is smaller than the dots; the error in the 
predicted point is a formal la error based on an extrapolation from the 
observations taken around the 2000 and 2005 perihelia. The vertical dashed 
lines in red show where the perihelia of 2000, 2005 and 201 1 fall. 

Figure 19. JPL torque model. Bold dashed segments depict constraints obtained 
from photometric fits, before, during and after the 2005 perihelion passage. 
Area of triangles are A= -74°, B= 68°, C= 245°. Area of rectangle D is -55°. 
Area A+B+D=49° is the estimate of the rotation runoff from the pre-perihelion, 
quiescent solution (Region B) at the epoch 31 days after perihelion, which 
corresponds to the Stardust-NExT flyby. 

Figure 20. Fit of the “Gauss” torque model to spin rate and rotational phase data 
spanning the 2005 perihelion and including the Deep Impact results. ( Top 
panel) The component torques (green and blue) and the net torque (red) 
scaled to unit moment of inertia as a function of time. Note that the net torque 
is predicted to maximize 32 days before perihelion passage (red dot-dash 
vertical line). The net torque is at first negative and then rises to a positive 
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maximum. The torque falls off becoming negative again before fading out. 

The bulk of the torque operates for ~ 300 days around perihelion. ( Bottom 
panel ) The modeled run of spin rate in the vicinity of perihelion. The observed 
spin rates are marked in blue and lie precisely on the model curve. In addition 
the rotational phase at each point is reproduced to within ~1°. The rate of 
change of the spin rate measured from the Deep Impact approach 
photometry is shown as an inclined dot-dash line. In this model the predicted 
slope near perihelion is not precisely satisfied with the predicted value ~ 10% 
less than the measured value. 

Figure 21 . The spin rates predicted by the “Gauss” model around the 201 1 
perihelion including the time of the Stardust-NExT encounter (Feb 14, 2011). 
The time of encounter is shown as a black dot-dash vertical line. The only 
measured spin rate in this figure is that at JD2455074 the other two are 
predictions. The spin rate at encounter is estimated as 213.47 °/d. See text for 
details of the prediction. The time of perihelion is denoted with a vertical 
dotted green line. 

Figure 22. Modeled evolution of the spin rate of 9P/Tempel 1 through three 
perihelion passages. The dashed blue line is the Gauss (see text), non- 
gravitational, torque model fit to the 2001-2004, Deep Impact, and post-2005 
spin rates and rotational phases. The predicted evolution of the spin rate 
through the encounter of the Stardust-NExT mission on Feb 14, 2011, is also 
shown. The black points are the spin rates deduced from the data in regions 
A,B, Deep Impact, and C that are listed in Table 2. The dashed line indicates 
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the slope determined from the Deep Impact approach photometry. The 
predicted spin rates at encounter on Feb 14, 201 1 (213.47 °/d) and near 
aphelion (212.52 °/d) following the 201 1 encounter are also shown. The 
vertical dashed lines in red show where the perihelia of 2000, 2005 and 201 1 
fall. 

Figure 23. Contour plot of the percentage of the 2005 Deep Impact imaging 
coverage that will be observable at Stardust-NExT encounter as a function of 
sub-solar longitude at the time-of-arrival and the spacecraft B-plane angle 
with respect to the sun-line. The blue area is <25%, the magenta area 
achieves 25-50%, the yellow area achieves 50-75%, and the grey area 
achieves 75-100%. The green dot denotes the arrival conditions for optimal 
viewing of the Deep Impact crater. The black dot and bar denotes the 
conditions for the nominal time-of-arrival on Feb 14, 201 1 20:42 given the 
predictions based on the observations analyzed here. Since it falls in the blue 
region of the plot the Level 1 NASA requirement of imaging at least 25% of 
the Deep impact coverage is unlikely to be attained without an adjustment of 
the time-of-arrival. A time-of-arrival adjustment of a delay of ~8 hr is required 
to move the black point and its estimated error into the magenta area and 
achieve the requirement. Such an adjustment was made by a trajectory 
correction maneuver in February, 2010. 

Figure 24. The anticipated view of the nucleus of 9P/Tempel 1 from the Stardust- 
NExT spacecraft at closest approach on a trajectory that has the sub- 
spacecraft latitude at -11°. Grey represents areas previously imaged by Deep 
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Impact while blue is terra incognita. Also shown is the position of the Deep 
Impact crater, the prime meridian, and the position of the South pole. The LH 
panel shows the view with the nominal time-of-arrival in January 2010; The 
RH panel shows the view with the maximum 8h delay in the time-of arrival. 

Figure A1 . A full frame MRI image showing the quadrant nomenclature. 

Figure A2: Illustration of the procedure used to create a uniform background in 
the MRI science images. 

Figure A3: Procedure used to correct comet flux measurements for the quadrant 
boundary effect (see text). The cross at the center of the aperture 
corresponds to the comet's centroid position. The green line represents the 
quadrant boundary. The correction consists of subtracting 1/3 of the flux from 
the rows at the boundary. The yellow shaded areas illustrate what is added 
back in. 

Figure B1. Dynamical periodograms calculated with the methods of Drahus and 
Waniak (2006) showing the variance, R, in the Deep Impact approach 
photometry for assumed values of the spin rate and the spin rate slope. The 
Top Panel shows a harmonics fit with 5 harmonics, and the Bottom Panel 
shows the fit using the DPDM with 60 bins and 5 covers. The grey scale to 
the left indicates the range of R. Note that in both cases the minimum 
variance falls well above dco/dt = 0. 


77 



Figure 1 


18 


I- 15 - 


12 


i i i i 

i 

1 I 1 l; 

i 

1 1 1 1 

i i 

i 

i 

• i 

' A i 

i 

• *i * ■ 

Vi.i 

i I'lM 

i 

i 

i 

i 

i 

■id 

Region A * ! f 

^ ■ i 

Region 0 ' 

f 

Region C 

i 

i 

i 

Region D 

(P re 2 000) j 

1 

l 

(Post2000) \ 

l 

L_ 

(P re 2 009) 

j (Post201 1) 

1 

1 


19 97 19 98 19 99 2000 2001 2 00 2 20 09 2004 2005 20 06 20 07 2008 2009 2010 2011 2012 2013 

Calendar Year 


78 




R1.1.0) 


Figure 3 


18 


17 - 


16 - 


15 - 


12 


. 1 1 

1 

1 1 1 1 

1 1 1 



(ft » 

. 1 = * 

• 

Region A * 
(P re 2 000) 

J 

■ > 

Region B 

(Post2000) * 

Region C 
(P re 2 009) 



t 

i— - 




18 ©7 18 8 8 18 8 8 20 0 0 2 001 2 00 2 2 00 3 2 00 4 2 00 5 2 00 6 2 00 7 2 00 8 2 00 8 2 0 1C 

Calendar Year 


Region D 
(Post2011) 


— 1 

2 0 11 2 0 12 2 0 18 


80 


Figu 



8 








Figure 6 



Calendar date 


83 


R(l.l.a) 


Y AU (ecliptic cooiri} 


Figure 7 



-3 -2 -1 0 1 2 3 4 5 


X AU (ecliptic coord) 


84 



Fraction of signal due to nucleus 


Figure 8 


1.0 
0.8 
0.6 
0.4 
0.2 
0.0 

-70 -60 -50 -40 -30 -20 -10 0 10 

Days before impact 



85 



Normalized Brightness in D 


Figure 9 



-60 -50 -40 -30 -20 -10 0 


Days before impact 


86 



Figure 10 





Cleaned-up Data (Brightness Units) 


Figure 11 



2001 JO 20015 20020 20025 20030 20035 20040 20045 20050 



♦ 

Region C 


20085 2007 20075 2008 20085 2009 20095 2010 20105 


Calendar Date 


88 





Figure 12 



89 


RMS of residuals 


Figure 13 



co - deg/d ay 


90 


Figure 14 




91 


















Figure 16 


60 • 
40 ■ 
20 • 
0 ■ 
-20 • 
-40 ■ 




' ; J&- : ■ 1 Jfc'- : a 


Region A 


c 

3 


<D 

£Z 

JZ 

n> 

CD 

03 

‘ 

03 

Q 

Q. 

13 

i 

C 

03 

_0J 

O 


0.0 0.2 0.4 0.6 0 JB 1.0 1.2 14 1.6 1 J8 2.0 




Phase 


93 


Figure 17 
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Figure 21 
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Figure 22 


Model Evolution of 9PfTempel 1 Spin Rate 
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Figure 23 
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Figure A2 
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Figure A3 


P.nt .lliove l>oun«Uiiy 



P.nt below bouml.iiy 



104 


dw/dt ( 1 O' 7 h'-) 


Figure B1 



00240 0.0245 0.0250 


(0 0 (h' 1 ) 


105 



